O 
O 



The U-shaped distribution of globular cluster specific frequencies 
in a biased globular cluster formation scenario 

Kenji Bekki 

School of Physics, University of New South Wales, Sydney NSW 2052, Australia 

Hideki Yahagi 



Department of Astronomy, University of Tokyo, 7-3-1 Hongo, Bunkyo ward, Tokyo 

113-0033, Japan 

and 



Duncan A. Forbes 



Centre for Astrophysics & Supercomputing, Swinburne University of Technology, 

Hawthorn, VIC 3122, Australia, 

oo 

\Q • ABSTRACT 
O 

o 



Using high-resolution numerical simulations, we investigate mass- and 
luminosity- normalized specific frequencies (Tn and Sn, respectively) of globular 



Q-f cluster systems (GCSs) in order to understand the origin of the observed U- 

shaped relation between Sn and U-band magnitude (My) of their host galaxies. 
We adopt a biased GC formation scenario in which GC formation is truncated 
in galaxy halos that are virialized at a later redshift, z trun . T$ is derived for 
galaxies with GCs today and converted into Sn for reasonable galaxy mass-to- 
^ ■ light-ratios (M/L). We find that Tn depends on halo mass (Mh) in the sense 

that Tn can be larger in more massive halos with Mh > 1O 9 M , if z tvun is as high 
as 15. We however find that the dependence is too weak to explain the observed 
Sn-Mv relation and the wide range of Sn in low-mass early-type galaxies with 
—20.5 < My < —16.0 mag for a reasonable constant M/L. The My-dependence 
of Sn for the low-mass galaxies can be well reproduced, if the mass-to-light-ratio 
Mh/ Ly oc M h a , where a is as steep as —1. Based on these results, we propose 
that the origin of the observed U-shaped Sn-Mv relation of GCSs can be under- 
stood in terms of the bimodality in the dependence of Mh/ Ly on Mh of their host 
galaxies. We also suggest that the observed large dispersion in Sn in low-mass 
galaxies is due partly to the large dispersion in T N . 



Subject headings: globular clusters: general - galaxies: star clusters - galaxies: 
evolution - galaxies: stellar content 
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1. Introduction 

Specific frequencies (Sn) of globular clusters (GCs) in galaxies are observed to be dif- 
ferent between galaxies with different Hubble morphological types and between those in 
different environments (e.g., Harris 1991). The Sn of globular cluster systems (GCSs) and 
their correlations with physical properties of their host galaxies have long been discussed 
in various different contexts of galaxy formation; biased GC formation in high-density en- 
vironments (e.g., West 1993), origin of field and cluster elliptical galaxies (e.g., Ashman & 
Zepf 1992; Forbes et al. 1997), and origin of dwarf elliptical (dE) galaxies (e.g., Miller et 
al. 1998). Recent numerical and theoretical studies have also discussed the origin of Sn of 
elliptical galaxies interacting with cluster environments (e.g., Bekki et al. 2003). 

Durrell et al. (1996) investigated the trend between Sn and M v (where M v is the 
U-band magnitude of a galaxy) for galaxies with a wide range of My. They found that 
Sn values were higher in the lowest luminosity dwarf ellipticals (dEs), and similar to the 
Sn values for giant ellipticals (see Fig. 9 in their paper). Sn values appeared to have a 
minimum at intermediate galaxy luminosities, but there were very few data points available 
to them. Forbes (2005) recently investigated the Sn values for ~100 early-type galaxies in 
the Virgo cluster, confirming the high Sn values for dwarfs and low Sn values at My ~ 
-19.5. Thus the S N -M V relation reveals a 'U'-shaped, or bimodal, distribution. Durrell et 
al. (1996) speculated that supernova-driven mass loss was responsible for the higher S N 
values in dwarfs. Forbes (2005) suggested that this, and the transition from hot to cold gas 
accretion (Dekel & Birnboim 2005) played a role in determining the number of GCs per unit 
galaxy starlight. He further noted that this increase in Sn closely followed the dependence 
of galaxy M/L with mass. However, this observed Sn bimodality has not yet been explored 
by theoretical/numerical studies of GC formation. 

The purpose of this Letter is to demonstrate, for the first time, whether and how the 
observed U-shaped Sn-Mv relation is reproduced in a "biased GC formation scenario" in 
which GCs can be formed only in galaxy subhalos with high primordial densities and thus 
virialized earlier (i.e., further GC formation is truncated after a redshift of z trmi ). This biased 
GC formation scenario has been discussed in the context of GCS properties (e.g., West 1993; 
Beasley et al. 2002; Santos 2003; Bekki 2005; Yahagi & Bekki 2005; Moore et al. 2005), 
but not discussed in the context of the U-shaped Sn distribution. Throughout this paper, 
we adopt the standard definition of S N = N GC x 10-°- 4 ( Mv+15 ) and T N = T N GC /M, where 
A^gc is the total number of GCs in a GCS and T is a normalization factor such that Tn has 
reasonable values (an order of 1 — 10). 
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2. The observed U-shaped relation and its fitting function 

Figure 1 shows the observed U-shaped S^-My relation from Harris (1991) for bright 
galaxies and from Durrell et al. (1996) for dwarf ones. It is clear from this figure that the 
galaxy luminosity dependence of is different between galaxies below and above a threshold 
luminosity, M Vit h- Although the dependence is less clear for galaxies with My < M Vjt h in 
Figure 1, Rhode et al. (2005) found a strong trend of increasing T N with increasing galaxy 
luminosity for elliptical galaxies with mass > 10 n M Q . 

The observed S^-My and T N -M relations may provide vital clues to galaxy and GC 
formation (e.g., Harris 1991; Forbes 2005). Accordingly comparison between observed and 
simulated relations may well enable us to derive some physical meaning for the relations. It 
is therefore important to quantify the relations by using an admittedly simple mathematical 
expression. We here propose the following form for the S^-My relation; 

S N (x) = A 1 x W K ^ + (S N ,th - A ± ) x 10^, (1) 

where x = Mv ^^' th and S Njth is S N at My = M Vit h- Parameter values of A 1: K 1: K 2 , M Vit h, 
and SV^th can be determined by fitting to observations. An example of the model fit with 
My :th = —19.5 mag, SN,th = 1-0, A\ = 0.5, K\ = —6, and K 2 = 4 is shown in Figure 1. 

This functional form in equation (1) can be equivalent to B\ x Ly Cl + B 2 x Ly ' 2 , where 
Ly is the U-band luminosity and other parameter values (e.g., B\) can be derived from 
values of the parameters (e.g., Ai) in equation (1). Therefore, the functional form has some 
advantages in discussing both theoretically and observationally the origin of 5n (and Tn) in 
terms of mass and luminosity growth processes of galaxy assembly. The observed 7 of 1 — 2 
in iVcc oc Ly 1 for bright ellipticals (Zepf & Ashman 1993) suggests that C 2 can range from 
to 1. 



3. Numerical simulations of T N and 
3.1. The model 

Since methods and techniques of numerical simulations of GC formation at high redshifts 
are given in detail in our previous paper (Yahagi & Bekki 2005, YB), we only briefly describe 
them in this paper. We simulate the large scale structure of GCs in a ACDM Universe 
with Vt = 0.3, A = 0.7, H = 70 km s -1 Mpc -1 , and <r 8 = 0.9 by using the Adaptive 
Mesh Refinement N— body code developed by Yahagi (2005) and Yahagi & Yoshii (2001). 
We use 512 3 collisionless dark matter (DM) particles in a simulation with the box size of 
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8.26/i -1 Mpc and a particle mass of 5 x 1O 5 M . We start simulations at z — 41 and follow 
them until z = in order to investigate the physical properties of old GCs outside and 
inside virialized dark matter halos at z — 0. We used the COSMICS (Cosmological Initial 
Conditions and Microwave Anisotropy Codes), which is a package of fortran programs for 
generating Gaussian random initial conditions for nonlinear structure formation simulations 
(Bertschinger 1995, 2001). 

We assume that (1) old, metal-poor globular cluster (MPC) formation is truncated for 
z < z trmi , and (2) initial radial profiles of GCSs in virialized subhalos at z — z trnn are the 
same as those of the "NFW" profiles (Navarro, Frenk, & White 1996). We follow dynamical 
evolution of GCs from z = z trm i to and thereby investigate mass of the virialized halos at 
z = (M h ), the total numbers of GC particles within the halos (N GC ), and the normalized 
number of GCs per halo mass (T N = T N GC /M h ). Given the fact that there are currently 
no observational constraints on initial GCS density profiles at z — z trun , we consider that 
adopting the NFW profiles is a reasonable first step for better understanding GCS properties 
in the CDM model. 

We convert the simulated T N -M h relations into the S N -My ones by using the following 
formula: 



where Sq is a constant and the normalization factor T for each model is determined such 
that number of GCs per unit mass in a simulation box is the same between models with 
different z trun . We then compare the simulated relations to the observed data of Forbes 
(2005) in which the of blue GCs was derived based on new data from the ACS Virgo 
Cluster Survey (Peng et al. 2005). 

We investigate two M/L models, i.e. variable and constant M/L. For the variable M/L 
model, we assume that M h /Ly oc M h a for M h < M hjth , where M hjth is set to be lO n M , and 
M h /L v oc Mh' 3 for M h > M h)th . We here focus mainly on low-luminosity galaxies for which 
observational data of Sn is available (Forbes 2005). We accordingly chose the value of (3 
consistent with observations based on galaxy luminosity functions (i.e. (3 = 0.33; Marinoni 
& Hudson 2002) and derive the best value of a for low-luminosity galaxies. Although we 
investigate models with a wide range of M^/Ly (i.e. 5 — 100) for the constant M/L model, 
we only show the results of the model with M^/Ly = 10. It should be stressed here that 
all of the present results are galaxies containing GCs at z = 0; GC-less galaxies at z = 
(typically those with My > —15) are not included in the present study. 

Although Forbes (2005) speculated about a possible M-dependent Mj L explaining the 
observed S^-My relation, his adopted assumptions (e.g., constant Tn) are not based on 
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detailed theoretical/numerical studies. Since the present study derives 5n from the simulated 
Tn, we can discuss more self-consistently the origin of the S^-My relation on firm physical 
basis. We demonstrate that the simulated T^-M^ dependences are too weak to explain the 
wide range of 5n for a constant M/L: The difference in GC number per unit galaxy mass 
between galaxies alone can not explain observations. 

3.2. Results 

Figure 2 shows the distributions of the simulated GCSs at z = in the T^-M^ plane for 
the models with z trun = 6 and 15. Below a critical halo mass of ~ 2 x 1O 9 M , Tn is higher for 
halos with smaller in the model with z trnn = 15. Above the critical mass, Tn is higher for 
halos with larger M n in this model and the slope of the Mh-dependence of Tn is significantly 
shallower than in halos below M n; th- This clear bimodality in the T-^-M^ relation is not seen 
in the model with z trmi = 6, which shows a monotonically increasing trend of T N with M h 
for 1O 8 M < M h < 3 x 10 10 M Q . The dispersion in T N appears to be significantly larger in 
halos with smaller M n for both models, which implies that the observed large dispersion in 
Sn in low-mass galaxies can be understood in terms of the large dispersion in Tn. 

For the model with z trun = 15, GC formation is assumed to occur only in the subgalactic 
halos formed from rare high-density peaks in the primordial matter distribution and thus 
virialized before z = 15. Accordingly, the bimodality can be understood in terms of the 
following two competing effects of halo merging with, or without, GCs. Low- mass subhalos 
with initially high Tn (virialized at z > 15) grow after their formation by merging/accretion 
of subhalos most of which are virialized later than z — 15 and thus are those without GCs 
(referred to as "GC-less" halos). Consequently, their Tn values can decrease as their masses 
become larger during hierarchical growth. This is the first effect of halo merging. 

According to the theory of biased galaxy formation (e.g., Kaiser 1984), the primordial 
density field of a more massive halo today is more likely to contain a larger fraction of high- 
density peaks. Accordingly, a larger halo at z = can be formed from a larger fraction of 
subhalos that are virialized before z = 15 and which contain GCs ("GC" halos). Therefore, 
a larger halo is more likely to be formed from merging of a larger fraction of GC halos and 
thus merging can increase T N for larger halos. This is the second effect of halo merging. 

As a result of these two competing effects, we can understand the bimodal T N -M h 
relation, i.e. below M^th the first effect is stronger, while above the second effect is stronger. 
The dramatic effect of the first GC-less halo merging is not expected for the model with 
-^trun = 6, because most of subgalactic halos that are building blocks of larger halos at z = 
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have been virialized before z = 6 and contain GCs. Therefore the monotonically increasing 
TN-Mh relation is due to the second effect in this model. It should be stressed here that 
irrespective of z tiun , T N is more likely to be larger for halos with M h > 1O 9 M . 

Figure 3 shows the comparison between the observations (Forbes 2005) and the best 
variable M/L model with aw-1, where a is the slope in the relation of M h /Ly oc M h Q . 
It should be noted here that the observed Sn of blue GCs in ellipticals brighter than My ~ 
—20.5 mag are seriously underestimated due to the limited ACS field-of-view (Forbes 2005). 
Thus the comparison makes sense only for lower luminosity galaxies. As shown in Figure 3, 
both models with z trnn = 6 and z tIUI1 = 15 can reproduce well the observed trend of increasing 
Sn with decreasing luminosity for low- mass galaxies with My > M^th- 

Figure 4 shows that the simulated S-^-My relations do not match well with the observa- 
tions for the constant M/L models. We also find that constant M^/Ly models, with a wide 
range of M^/Ly values (i.e. 5 — 100), do not match with the observations either. These re- 
sults, combined with those in Figure 3, imply that the origin of the S^-M^ relation is closely 
associated with the dependence of M h / Ly on M h , as suggested by Forbes (2005). Figure 4 
shows that Sn can be only slightly higher in low-luminosity galaxies with My > —16.5 mag 
as a result of biased GC formation in the model with z tT nn = 15. However the simulated 
dependence is not strong enough to explain quantitatively the observation for the constant 
Mh/Ly model. Since observations for galaxies brighter than My = —16 mag should be com- 
pared with simulations, the results for Mh > 2 x 1O 9 M in Figure 2 are relevant to Figure 
4. The derived S^-My relations in Figure 4 are thus essentially the same as T^-M^ ones for 
Mh > 2 x 1O 9 M in Figure 2. If we adopt a higher M/L, the blue/red lines are shifted to 
the left in Figure 4: A model with M/L = 25 results in a 1 mag shift to lower luminosities. 

4. Discussion and conclusions 

West (1993) first suggested that the observed dependence of Sn on galaxy type and 
environment can be understood in terms of GC formation from the statistical enhancement 
of rare high-density peaks in the primordial matter distribution (i.e., biased GC formation). 
The present numerical simulations have shown that biased GC formation in subgalactic halos 
virialized early (z ~ 15), can be imprinted on the observed T N — M h relation. Semi-analytic 
galaxy formation models based on a hierarchical clustering scenario have also shown that 
models with a truncation of GC formation at z ~ 5 can better reproduce the observed color 
bimodality of GCs in early-type galaxies (Beasley et al. 2002). What physical mechanisms 
are responsible for the truncation (or the severe suppression) of GC formation in subgalactic 
halos virialized at later times? 
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Santos (2003) first proposed that cosmic reionization at high z can suppress GC for- 
mation and thus be an important factor for better understanding the origin of the physical 
properties of GCSs. Recent numerical simulations on GCs and old stellar halos in galaxies 
have confirmed that structural and kinematical properties of GCSs carry fossil information 
about the epoch of the reionization -2 re ion that can suppress GC formation (Bekki 2005; Bekki 
& Chiba 2005; Moore et al. 2005; Rhode et al. 2005). Therefore, if z trU n = -Zrcion, the present 
results imply that the T N — M h relation can provide some constraints on z rcion . Furthermore, 
we suggest that these relations can be different in different environments (or cosmic volumes 
investigated by observations), if reionization is spatially inhomogeneous (Benson et al. 2001; 
Wyithe & Loeb 2004). 

We suggest that it is worthwhile for future observations to investigate the T N -M h relation 
for a statistically significant number of galaxies, although it is a formidable task to derive 
total galactic masses including dark matter halos. Such observational studies can test biased 
GC formation at high redshifts by comparing the results with the corresponding simulations. 
The predicted Tn values change only by a factor of ~ 5 for a change in Mh by four orders 
of magnitude. The observed difference in Sn for a My range of ~ 10 mag (corresonding 
roughly to a luminosity difference of factor ~ 10 4 ) is more than a factor of ~ 30 (see Figure 
1). Therefore the observed bimodal S^-My relation can not be explained well by the present 
models without resorting a mass-dependent Mj L (for which there is some independent evi- 
dence). We accordingly suggest that the origin of the observed U-shaped S^-My relation of 
GCSs can be associated, at least partly, with the bimodality in the mass-dependence of the 
Mj L ratio of their host galaxies. 
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to improve the present paper. The numerical simulations reported here were carried out on 
Fujitsu-made vector parallel processors VPP5000 kindly made available by the Astronomical 
Data Analysis Center (AD AC) at National Astronomical Observatory of Japan (NAOJ) for 
our research project why36b. K.B. and D.A.F. acknowledge the financial support of the 
Australian Research Council throughout the course of this work. H.Y. acknowledges the 
support of the research fellowships of the Japan Society for the Promotion of Science for 
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Fig. 1. — Observational data for GCSs in galaxies in the S^-My plane. Observational 
data are from Harris (1991) for bright galaxies and from Durrell et al. (1996) for dwarf 
populations. The solid line represents an empirical fit to the observed bimodality in the 
relation between Sn and My (the details of the functional form are given in the main text). 
The dotted line marks the critical magnitude (M Vi th) : The dependence on My is different 
between low and high luminosity galaxies. 
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Fig. 2. — Dependences of T N (=T N GC /M h ) on log 10 M h for the models with z trmi = 6 (upper) 
and z t run = 15 (lower). Small dots represent simulated GCSs within galaxy-scale halos of 
masses M h and "GC particle" numbers of N GC - Big dots with error bars represent the mean 
values of T N in each mass bin. The error bar in each mass bin is estimated as T N /-y/2(A^ — 1), 
where iVj is the total number of galaxy-scale halos in the i-th mass bin. In order to compare 
these results with observations, the value of the normalization factor (T ) is determined for 
each model. See text for details. 
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Fig. 3. — Observational data for ~100 GCSs in early-type galaxies in the S^-My plane. 
Here Sn values are shown for blue GCs only from Forbes (2005) and the large triangles show 
the mean value of Sn in each magnitude bin. For comparison, the variable M/L models 
with Mh/Ly oc Mh -1 ' are shown for z trun = 6 (solid) and z t run = 15 (dotted). Since Sn 
of brighter galaxies with My < —20.5 mag are significantly underestimated in observations, 
the lower limits of Sn values are indicated by green arrows for these galaxies. 
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Fig. 4. — The same as Figure 3 but for the models with constant M^/Ly (=10). The shapes 
of the derived Sn — My relation can be slightly different from those of the — M n one 
shown in Figure 2, though overall shapes are quite similar between the two relations for 
Mh > 8.5 x 1O 8 M (corresponding to My < —15.0 mag). The origin of the slight differences 
is associated with the difference in the binning process of simulation data between Figures 
2 and 4. 



